Setting up

## Loading required package: maptools
## Loading required package: sp
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
## Loading required package: RColorBrewer
## Loading required package: MASS
## Loading required package: rgeos
## rgeos version: 0.5-9, (SVN revision 684)
##  GEOS runtime version: 3.9.1-CAPI-1.14.2 
##  Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
##  GEOS using OverlayNG
##  Linking to sp version: 1.5-0 
##  Polygon checking: TRUE
## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
## 
## Attaching package: 'raster'
## The following object is masked from 'package:MASS':
## 
##     select
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/user/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/user/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 2.4-0
## 
## Attaching package: 'spatstat.geom'
## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift
## The following object is masked from 'package:MASS':
## 
##     area
## Loading required package: spatstat.random
## spatstat.random 2.2-0
## Loading required package: spatstat.core
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:raster':
## 
##     getData
## Loading required package: rpart
## spatstat.core 2.4-4
## Loading required package: spatstat.linnet
## spatstat.linnet 2.3-2
## 
## spatstat 2.3-4       (nickname: 'Watch this space') 
## For an introduction to spatstat, type 'beginner'

Loading data

Polygons and lines

planning_area_sf <- st_read("data/Planning_Area_Census2010.kml")
## Reading layer `Planning_Area_Census2010' from data source 
##   `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\Planning_Area_Census2010.kml' 
##   using driver `KML'
## Simple feature collection with 55 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XYZ
## Bounding box:  xmin: 103.6054 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
class(planning_area_sf)
## [1] "sf"         "data.frame"

Points

Next we read in point data. The first set of point data is the Park Connector (PCN) access points themselves. Next, we have data on various public amenities like supermarkets, hawker centres, pharmacies as well as transportation infrastructure like MRT exits.

access_points_sf <- st_read("data/PCN_Access_Points.kml")
## Reading layer `PCN_Access_Points' from data source 
##   `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\PCN_Access_Points.kml' 
##   using driver `KML'
## Simple feature collection with 306 features and 2 fields
## Geometry type: POINT
## Dimension:     XYZ
## Bounding box:  xmin: 103.6933 ymin: 1.272978 xmax: 104.0048 ymax: 1.46272
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
supermarkets_sf <- st_read("data/supermarkets.kml")
## Reading layer `SUPERMARKETS' from data source 
##   `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\supermarkets.kml' 
##   using driver `KML'
## Simple feature collection with 742 features and 2 fields
## Geometry type: POINT
## Dimension:     XYZ
## Bounding box:  xmin: 103.6258 ymin: 1.24715 xmax: 104.0042 ymax: 1.462167
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
hawkercentres_sf <- st_read("data/hawkercentre.kml")
## Reading layer `HAWKERCENTRE' from data source 
##   `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\hawkercentre.kml' 
##   using driver `KML'
## Simple feature collection with 125 features and 2 fields
## Geometry type: POINT
## Dimension:     XYZ
## Bounding box:  xmin: 103.6974 ymin: 1.272716 xmax: 103.9882 ymax: 1.449217
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
mrt_exits_sf <- st_read("data/lta-mrt-station-exit-kml.kml")
## Reading layer `MRT_EXITS' from data source 
##   `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\lta-mrt-station-exit-kml.kml' 
##   using driver `KML'
## Simple feature collection with 474 features and 2 fields
## Geometry type: POINT
## Dimension:     XYZ
## Bounding box:  xmin: 103.6368 ymin: 1.264972 xmax: 103.9893 ymax: 1.449157
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
pharmacy_sf <- st_read("data/registered_pharmacy.kml")
## Reading layer `REGISTERED_PHARMACY' from data source 
##   `C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data\registered_pharmacy.kml' 
##   using driver `KML'
## Simple feature collection with 269 features and 2 fields
## Geometry type: POINT
## Dimension:     XYZ
## Bounding box:  xmin: 103.699 ymin: 1.2563 xmax: 104.0025 ymax: 1.448227
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
class(access_points_sf)
## [1] "sf"         "data.frame"
class(supermarkets_sf)
## [1] "sf"         "data.frame"
class(hawkercentres_sf)
## [1] "sf"         "data.frame"
class(mrt_exits_sf)
## [1] "sf"         "data.frame"
class(pharmacy_sf)
## [1] "sf"         "data.frame"

Raster

Next, we read in raster data, which consists of only population density data.

GDALinfo("data/popmap15adj.tif")
## Warning in GDALinfo("data/popmap15adj.tif"): statistics not supported by this
## driver
## rows        368 
## columns     542 
## bands       1 
## lower left origin.x        103.6383 
## lower left origin.y        1.164905 
## res.x       0.0008333 
## res.y       0.0008333 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=longlat +datum=WGS84 +no_defs 
## file        data/popmap15adj.tif 
## apparent band summary:
##    GDType hasNoDataValue   NoDataValue blockSize1 blockSize2
## 1 Float32           TRUE -3.402823e+38        128        128
## apparent band statistics:
##          Bmin       Bmax Bmean Bsd
## 1 -4294967295 4294967295    NA  NA
## Metadata:
## AREA_OR_POINT=Area
pop_density <- raster("data/popmap15adj.tif")
pop_density
## class      : RasterLayer 
## dimensions : 368, 542, 199456  (nrow, ncol, ncell)
## resolution : 0.0008333, 0.0008333  (x, y)
## extent     : 103.6383, 104.09, 1.164905, 1.471559  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : popmap15adj.tif 
## names      : popmap15adj

Exploratory Data Analysis

Park Connectors

First, we plot the PCN access points. On first glance, it is apparent that these points cannot be randomly distributed, with noticeable clusters in the Choa Chu Kang/Bukit Batok area as well as the Woodlands/Sembawang area. The points also seem to occur as a “string” of points, which obviously stands to reason since these are access points to Park Connectors which are essentially lines.

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_borders() + tmap_options(check.and.fix = TRUE) + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "red")+ tm_layout(title="PCN Access Points in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid

Plot various amenities over the PCN access points

Now, we will display the geographical distribution of the various amenities (supermarkets, hawker centres, etc.) on top of the locations of the PCN acces points. ### Supermarkets We plot the locations of supermarkets overlaid with the access points on a map of Singapore. We see that most access points are located very close to supermarkets.

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_borders() + tmap_options(check.and.fix = TRUE) + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "red") + tm_shape(supermarkets_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "skyblue") + tm_layout(title="PCN Access Points in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid

Hawker centres

For hawker centres, the correlation is not as obvious.

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_borders() + tmap_options(check.and.fix = TRUE) + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "red") + tm_shape(hawkercentres_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "green") + tm_layout(title="Cycling paths in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid

MRT exits

For MRT exits, we notice that the Downtown Core area has a high concentration of exits. This could be because many of the stations in that area have many exits, sometimes as many as 7 or 8. This makes it seem more numerous. We can try getting a dataset of just MRT stations themselves to help us narrow down our investigation. Unlike the MRT exits, we see that there are very few PCN access points in the Downtown Core area. This could be because of the lower priority assigned to cycling paths over mass rapid transit in a space-scarce area like the Downtown Core.

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_borders() + tmap_options(check.and.fix = TRUE) + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "red") + tm_shape(mrt_exits_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "orange") + tm_layout(title="Cycling paths in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid

Pharmacies

For pharmacies, there doesn’t seem to be a clear pattern.

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_borders() + tmap_options(check.and.fix = TRUE) + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "red") + tm_shape(pharmacy_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "purple") + tm_layout(title="Cycling paths in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid

Population density

Interestingly, areas with the highest population density (Serangoon, Hougang and the areas surrounding Paya Lebar) seem to be relatively underserved by acces points as compared to Western regions like Bukit Batok which has a much lower population density.

Standing out as well is the central part of Singapore, which has virtually zero population density. This region encompasses the MacRitchie Reservoir area and naturally is not home to any human population density. It might be useful to exclude this region in later analysis to prevent confounding.

tmap_mode("view")
## tmap mode set to interactive viewing
# tm_shape(pop_density) + tm_raster() + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "green") + tm_layout(title="Cycling paths in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))

tm_shape(pop_density) + tm_raster() + tm_shape(access_points_sf) + tm_symbols(size = 0.1, alpha = 0.5, col = "green") + tm_shape(planning_area_sf) + tm_text("Name", size=0.5) + tm_layout(title="Cycling paths in Singapore", title.size = 0.65, title.position = c(0.02, "bottom"))
## Warning: The shape planning_area_sf is invalid. See sf::st_is_valid

Extracting planning areas to observe

training_sf <- planning_area_sf[planning_area_sf$Name %in% c("YISHUN","JURONG EAST", "JURONG WEST", "TAMPINES", "PASIR RIS", "BEDOK"),]
training_sf
## Simple feature collection with 6 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XYZ
## Bounding box:  xmin: 103.6749 ymin: 1.296503 xmax: 103.9848 ymax: 1.445616
## z_range:       zmin: 0 zmax: 0
## Geodetic CRS:  WGS 84
##           Name
## 1    PASIR RIS
## 8  JURONG EAST
## 14 JURONG WEST
## 23      YISHUN
## 35    TAMPINES
## 52       BEDOK
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             Description
## 1       <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>PASIR RIS</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>PASIR RIS</td> </tr> <tr> <td>Planning Area Code</td> <td>PR</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>EAST REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>ER</td> </tr> <tr> <td>INC_CRC</td> <td>77F6CE717F532D1A</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:19 PM</td> </tr> <tr> <td>X_ADDR</td> <td>40795.475</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>40066.2774</td> </tr> <tr> <td>SHAPE_Length</td> <td>22444.999809</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>15810712.053754</td> </tr> </table> </td> </tr> </table> </body> </html>
## 8   <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>JURONG EAST</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>JURONG EAST</td> </tr> <tr> <td>Planning Area Code</td> <td>JE</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>WEST REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>WR</td> </tr> <tr> <td>INC_CRC</td> <td>E0C3B6AB25F8AFFF</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:19 PM</td> </tr> <tr> <td>X_ADDR</td> <td>17035.1876</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>33633.557</td> </tr> <tr> <td>SHAPE_Length</td> <td>27342.347259</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>17857762.682791</td> </tr> </table> </td> </tr> </table> </body> </html>
## 14 <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>JURONG WEST</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>JURONG WEST</td> </tr> <tr> <td>Planning Area Code</td> <td>JW</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>WEST REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>WR</td> </tr> <tr> <td>INC_CRC</td> <td>3302EAE21AD1E998</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:20 PM</td> </tr> <tr> <td>X_ADDR</td> <td>13705.4978</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>35973.1233</td> </tr> <tr> <td>SHAPE_Length</td> <td>19356.680643</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>14688294.312111</td> </tr> </table> </td> </tr> </table> </body> </html>
## 23           <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>YISHUN</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>YISHUN</td> </tr> <tr> <td>Planning Area Code</td> <td>YS</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>NORTH REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>NR</td> </tr> <tr> <td>INC_CRC</td> <td>8D4108597B11C0DC</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:19 PM</td> </tr> <tr> <td>X_ADDR</td> <td>28472.3369</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>44137.8596</td> </tr> <tr> <td>SHAPE_Length</td> <td>20778.98769</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>21578985.748052</td> </tr> </table> </td> </tr> </table> </body> </html>
## 35       <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>TAMPINES</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>TAMPINES</td> </tr> <tr> <td>Planning Area Code</td> <td>TM</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>EAST REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>ER</td> </tr> <tr> <td>INC_CRC</td> <td>C5A6D638AA4E65FE</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:20 PM</td> </tr> <tr> <td>X_ADDR</td> <td>41454.9098</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>36208.5982</td> </tr> <tr> <td>SHAPE_Length</td> <td>22402.760996</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>20913646.028123</td> </tr> </table> </td> </tr> </table> </body> </html>
## 52             <html xmlns:fo="http://www.w3.org/1999/XSL/Format" xmlns:msxsl="urn:schemas-microsoft-com:xslt"> <head> <META http-equiv="Content-Type" content="text/html"> <meta http-equiv="content-type" content="text/html; charset=UTF-8"> </head> <body style="margin:0px 0px 0px 0px;overflow:auto;background:#FFFFFF;"> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-collapse:collapse;padding:3px 3px 3px 3px"> <tr style="text-align:center;font-weight:bold;background:#9CBCE2"> <td>BEDOK</td> </tr> <tr> <td> <table style="font-family:Arial,Verdana,Times;font-size:12px;text-align:left;width:100%;border-spacing:0px; padding:3px 3px 3px 3px"> <tr> <td>SHAPE</td> <td>Polygon</td> </tr> <tr bgcolor="#D4E4F3"> <td>Planning Area Name</td> <td>BEDOK</td> </tr> <tr> <td>Planning Area Code</td> <td>BD</td> </tr> <tr bgcolor="#D4E4F3"> <td>Central Area Indicator</td> <td>N</td> </tr> <tr> <td>Region Name</td> <td>EAST REGION</td> </tr> <tr bgcolor="#D4E4F3"> <td>Region Code</td> <td>ER</td> </tr> <tr> <td>INC_CRC</td> <td>DBC1F8A651648374</td> </tr> <tr bgcolor="#D4E4F3"> <td>FMEL_UPD_D</td> <td>4/14/2014 2:43:21 PM</td> </tr> <tr> <td>X_ADDR</td> <td>38582.8088</td> </tr> <tr bgcolor="#D4E4F3"> <td>Y_ADDR</td> <td>34032.1699</td> </tr> <tr> <td>SHAPE_Length</td> <td>21853.967373</td> </tr> <tr bgcolor="#D4E4F3"> <td>SHAPE_Area</td> <td>21731365.846805</td> </tr> </table> </td> </tr> </table> </body> </html>
##                          geometry
## 1  MULTIPOLYGON Z (((103.9532 ...
## 8  MULTIPOLYGON Z (((103.7119 ...
## 14 MULTIPOLYGON Z (((103.7281 ...
## 23 MULTIPOLYGON Z (((103.8468 ...
## 35 MULTIPOLYGON Z (((103.9643 ...
## 52 MULTIPOLYGON Z (((103.9636 ...

Calculating global density

sg_area <- 728.6 #km2
supermarkets_gd <- 742 / sg_area
hawkercentres_gd <- 125 / sg_area
mrt_exits_gd <- 474 / sg_area
pharmacy_gd <- 269 / sg_area

Quadrat density

sg <- readOGR("./data", "Region_Census2010")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\user\Desktop\SL Internship\BT4015\bt4012-geo-proj\data", layer: "Region_Census2010"
## with 5 features
## It has 8 fields
sg
## class       : SpatialPolygonsDataFrame 
## features    : 5 
## extent      : 2637.869, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs 
## variables   : 8
## names       : OBJECTID,       REGION_N,          INC_CRC, FMEL_UPD_D,    X_ADDR,     Y_ADDR,    SHAPE_Leng,    SHAPE_Area 
## min values  :        1, CENTRAL REGION, 51D899B637F87BAF, 2014/05/06, 13007.354, 31931.7133, 60174.6538079, 112912795.479 
## max values  :        5,    WEST REGION, C05F166F6716C12E, 2014/05/06,  42244.19, 44177.9547,  238296.73891, 251603415.909
sg_window <- as.owin(sg)
# convert crs of access_points_sf to crs of sg
access_points_q <- st_transform(access_points_sf, crs(sg))
# convert sf to ppp
access_points_ppp <- as.ppp(st_coordinates(access_points_q), W=sg_window)
## Warning: data contain duplicated points

Try plotting quadrats of various sizes.

Q6 <- quadratcount(access_points_ppp, nx=6, ny=6)
plot(Q6)

Q10 <- quadratcount(access_points_ppp, nx=10, ny=10)
Q15 <- quadratcount(access_points_ppp, nx=15, ny=15)
Q20 <- quadratcount(access_points_ppp, nx=20, ny=20)
Q40 <- quadratcount(access_points_ppp, nx=40, ny=40)
plot(Q40)

From the quadrat density graph, we again see that most access points are concentrated in the North and West.

# compute density within each quadrat
Q_density <- intensity(Q40)

# plot the density raster
plot(intensity(Q10, image=TRUE), main=NULL, las=1)
# plot access points
plot(access_points_ppp, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)

# plot the density raster
plot(intensity(Q20, image=TRUE), main=NULL, las=1)
# plot access points
plot(access_points_ppp, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)

# plot the density raster
plot(intensity(Q40, image=TRUE), main=NULL, las=1)
# plot access points
plot(access_points_ppp, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)

Quadrat density based on population density

pop_density_q <- projectRaster(pop_density, crs = crs(sg))

brks <- c(0, 50, 100, 150, 200, 250, 300)
zcut <- cut(pop_density_q, breaks = brks)
pop_tess <- tess(image = zcut)

Q2 <- quadratcount(access_points_ppp, tess = pop_tess)
## Warning in quadratcount.ppp(access_points_ppp, tess = pop_tess): Tessellation
## does not contain all the points of X
plot(Q2, main="", las=1)

# compute density within each quadrat
Q2_density <- intensity(Q2)
Q2_density
## tile
##            1            2            3            4            5            6 
## 3.029003e-07 1.125329e-06 1.760026e-06 5.016187e-07 4.759687e-07 5.626379e-07
cl <-  interp.colours(c("blue", "green", "yellow", "orange", "red", "darkred"), 6)
# plot the density raster
plot(intensity(Q2, image=TRUE), main=NULL, las=1, col=cl)
# plot access points
plot(access_points_ppp, pch=20, cex=0.6, col=rgb(0,0,0,.5), add=TRUE)

# Monte Carlo Q20 (Shu Ling)

ann_observed <- mean(nndist(access_points_ppp, k=1))
ann_observed
## [1] 300.9152

Population density

n <- 1000 # number of simulations
ann_pop <- vector(length = n) # empty object to be used to store simulated ANN values
pop_density_im <- as.im(pop_density_q) # convert population density to im object

for (i in 1:n) {
  random_point <- rpoint(n=access_points_ppp$n, f=pop_density_im)
  ann_pop[i] <- mean(nndist(random_point, k=1))
}

Window(random_point) <- sg_window
plot(random_point, pch=16, main=NULL, cols=rgb(0,0,0,0.5))

hist(ann_pop, main=NULL, las=1, breaks=40, col="bisque", xlim=range(ann_observed, ann_pop))
abline(v=ann_observed, col="blue")

## MRT exits

# convert crs of access_points_sf to crs of sg
mrt_exits_q <- st_transform(mrt_exits_sf, crs(sg))
# convert sf to ppp
mrt_exits_ppp <- as.ppp(st_coordinates(mrt_exits_q), W=sg_window)

MRT_Q20 <- quadratcount(mrt_exits_ppp, nx=20, ny=20)
plot(MRT_Q20)

MRT_density20 <- intensity(MRT_Q20, image=TRUE)
plot(MRT_density20)

class(MRT_density20)
## [1] "im"
n <- 1000 # number of simulations
ann_mrt <- vector(length = n) # empty object to be used to store simulated ANN values

for (i in 1:n) {
  random_point <- rpoint(n=access_points_ppp$n, f=MRT_density20)
  ann_mrt[i] <- mean(nndist(random_point, k=1))
}

Window(random_point) <- sg_window
plot(random_point, pch=16, main=NULL, cols=rgb(0,0,0,0.5))

hist(ann_mrt, main=NULL, las=1, breaks=40, col="bisque", xlim=range(ann_observed, ann_mrt))
abline(v=ann_observed, col="blue")